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We consider a two-state model of non-Markovian stochastic resonance (SR) within the framework 
of the theory of renewal processes. Residence time intervals are assumed to be mutually independent 
and characterized by some arbitrary non- exponential residence time distributions which are mod- 
ulated in time by an externally applied signal. Making use of a stochastic path integral approach 
we obtain general integral equations governing the evolution of conditional probabilities in the pres- 
ence of an input signal. These novel equations generalize earlier integral renewal equations by Cox 
and others to the case of driving- induced non-stationarity. On the basis of these new equations a 
response theory of two state renewal processes is formulated beyond the linear response approxi- 
mation. Moreover, a general expression for the linear response function is derived. The connection 
of the developed approach with the phenomenological theory of linear response for manifest non- 
Markovian SR put forward in [ I. Goychuk and P. Hanggi, Phys. Rev. Lett. 91, 070601 (2003)] is 
clarified and its range of validity is scrutinized. The novel theory is then applied to SR in symmetric 
Oh! non-Markovian systems and to the class of single ion channels possessing a fractal kinetics. 
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I. INTRODUCTION 



> 

o 

0^ ' The concept of Stochastic Resonance (SR) has been originally put forward in order to explain the periodicity 
, of glacial recurrences on the Earth [J. It has gained, however, an immense popularity in the context of signal 
■ transduction in nonlinear stochastic systems in physics and biology |^ |^ |^ . Paradoxically enough, the detection 
— ' of beneficial input signals in the background stochastic fluctuations of a signal-transmitting physical system can be 
improved upon corrupting the information-carrying signal with input noise, or upon raising the level of intrinsic 
thermal noise. A first example of SR has been given for a continuous state bistable dynamics agitated by the thermal 
I . noise and periodically modulated by an external signal ilj . There exists a huge number of systems in physics, chemistry, 
T3 ■ and biology which do exhibit SR These range from the classical systems to the systems with distinct quantum 

Ch features P- 

^ Experimentally, SR has been demonstrated in various macroscopic systems, see, e.g., in the reviews 0iQ and the 

. . [ references therein. For a mesoscopic system containing a finite number of molecules SR has been first demonstrated 
^ ' experimentally in Q. The mesoscopic system in 6j consists of dynamically self-assembled alamethicin ion channels 
of variable size that are placed in a lipid membrane. Up to this date, however, there remains the challenge to 
demonstrate SR on the level of single stable molecules. Ion channels of biological membranes 0, present one 
' of the most appealing objects for such single-molecular studies. The invention of patch clamp technique by Neher 
and Sakmann in 1976 8] made such investigations possible. The single-molecular SR experiments which have been 
performed under the conditions of variable intrinsic thermal noise intensity 9], did not arrive at the convincing 
conclusions. A recent theoretical study jlOl | suggested a parameter regime where SR effect should indeed occur for 
a Shaker K+ channel under physiological conditions when external noise is added to the signal. This issue has 
further been examined theoretically in The present status calls for both theoretical, and experimental further 
investigations. Particularly, the presence of distinct memory effects in the dynamics of such single molecules as ion 
channels ^3 constitutes a major theoretical challenge. The non-Markovian features caused by these memory effects 
may be crucial for the occurrence of Stochastic Resonance on the level of single molecules. 

The gross features of the observed bistable dynamics can be captured by a two-state stochastic process x{t) that 
switches forth and back between two values Xi and X2 at random time points {ti}. Such a two-state random process 
can be directly extracted from filtered experimental data and then statistically analyzed. Basically, the process x{t) 
is characterized as follows: The sojourn in the state xi alternates randomly at ti into the sojourn in the state X2, 
then x{t) switches back to xi at time ti-i-i, and so on. If the sojourn time intervals Tj = tj+i — ti are independently 
distributed (a condition which we shall assume throughout the following), such two-state renewal processes are fully 
specified by two residence time distributions (RTDs) tpi,2{T) jl^j- In the simplest case, which corresponds to the 
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dichotomic Markovian process (DMP), both RTDs are strictly exponential, i.e., iPi.2{t) = 1^1.2 exp(— z/i^2t), where z^i^2 
arc the transition rates which equal the inverse mean residence times (MRTs), which are given by 

poo 

(ti,2) := / rV-i.a^dr, (1) 
Jo 

with 1/1^2 — {ti^2)^^- The input signal f{t) causes the transition rates 1^1,2 to be time-dependent, i.e., ^1^2 — * '^i,2{t)- 
Moreover, the RTDs become functionals of the driving signal, 

^i^2{t't')^^Ji^2{t,t') = iyi,2{t)exp[- f iyi,2{T)dT]. (2) 



As a consequence, the time-dependent probabilities pi.2{t) of the states xi.2 obey the master equations 

P2{t) = Vl{t)pi{t) ~ V2{t)p2{t) (3) 

with the signal-dependent rates which under an adiabatic assumption obey the rate law jl4l |. 

^i,2(i) = v^il exp ( - [AC/1,2 T Axf{t)l2]lkBT) . (4) 



In 0), are the frequency prefactors, AC/1,2 are the heights of the activation barriers, Aa; :— X2 — xi > is the 
amplitude of fluctuations, ks is the Boltzmann constant and T is the temperature. For a weak periodic signal 

f{t)^focos{nt), (5) 

the use of Eqs. JSJ, Q allows one to calculate within linear response theory the asymptotic, long time response of 
the mean value {x{t)) = xipi{t) + X2P2{t) to f{t), i.e. 

{Sx{t)) ^ fo\x{fl)\cos{nt - (p{n)), ast->oo. (6) 

In JH)), is the linear response function in the frequency domain and ip{i^) denotes the phase shift. The spectral 

amphfication of signal, r] = |x(ri)|^, exhibits the effect of SR, i.e. a bell-shaped dependence versus increasing intrinsic 
thermal noise strength which is measured by the temperature T 0. 

The above outlined two-state Markovian theory has been put forward by McNamara and Wiesenfeld this 
approach has proven very useful over the years as a basic, prominent model for SR research 0- Remarkably enough, 
this simple model allows one to unify the various kinds of SR - periodic, aperiodic [l^ . and even non-stationary SR 
- within a unifying framework of information theory (lOj . 

Many observed bistable stochastic processes x{t) are, however, truly not Markovian, as can be deduced from the 
experimentally observed RTDs. As a matter of fact, any deviation of RTDs from the strictly exponential form indicates 
a deviation from the Markovian behavior [Hill. The profoundly non-Markovian case emerges when at least one of 
the RTDs possesses a large (diverging) variance var(ri,2) = /q°° T^^i^2iT)dT — (ti,2)^ 00. The stochastic dynamics 
of single molecules is especially interesting in this respect. For example, the RTDs of the conductance fluctuations 



example, tr 

in biological ion channels are in many cases not exponential [T^ I20I l2ll |. Usually, a sum of many exponentials, 
= IliLi ^i^i exp(— i^iT), X^ili Ci = 1 is needed to describe the experimental data 0. Moreover, in some cases 



V'(t) can well be described by a stretched exponential 1131 . o r by a power law ipir) oc \/{h + T)^ , /3 > |23,|21|. The 
power law is especially remarkable. For example, in Ref. [23 such a power law behavior has been found for the closed 
time RTD of a large conductance (BK) potassium channel with a power law exponent j3 « 2.24 yielding formally 
Yai{Tciosed) = oo. This in turn implies that such conductance fluctuations should exhibit a characteristic 1//" noise 
power spectrum S{f) jl^. Indeed, this is the case of BK ion channel [2^, as well as of some other ion channels [2^. 

What are the non-Markovian features of SR in similar systems? We address this question below using the 
just described non-Markovian generalization of McNamara- Wiesenfeld model characterized by some arbitrary non- 
exponential RTDs "01,2 (t) and the corresponding survival probabilities <I'i,2(''') = /°° i^i,2{T')dT' [l3- Similar models 
with alternating renewal processes have been used previously in the SR theory for some particular stochastic dynamics 
contracted to the two-state dynamics [l^ |2^ . Moreover, the class of colored noise driven Stochastic Resonance |23| 
is also intrinsically non-Markovian. All these prior studies have been restricted, however, to situations with finite 
memory effects on a finite time scale. A truly non-Markovian situation emerges when the memory effects extend 
practically to infinity, exhibiting a scale free, weak power law decay. A phenomenological linear response theory of 
such genuine non-Markovian SR (which does not presume a knowledge of the underlying microscopic dynamics) has 
been put forward recently in Ref. 28\ . The present work provides further details and, additionally, presents a more 
general framework for the non-Markovian SR theory which extends beyond the linear response description. 



3 



II. GENERAL THEORY 
A. Two-state renewal process 



To start, let us consider a two-state renewal process (TSRP) x{t) which takes initially, at time the value Xi, or 
the value X2 with the probability pi(io), or ^2(^0)7 correspondingly. At a random time point ti the process switches 
its current state with the probability one into another state and stays there until some next random time point t2. 
Then, the renewal process proceeds further in time in the same manner. The survival probability to remain in the 
state 1, or the state 2 for the time n = ti+j—ti is $i{Ti), or $2(7'i)> correspondingly. These two survival probabilities 
completely specify the considered TSRP T^. The functions ^i,2{t) must satisfy the following obvious restrictions: 
(i) < $i,2(t) < 1, (ii) ^i,2{t + Ar) < $i,2(t). At > (non-increasing function of time) (iii) $1,2(0) = 1; (iv) 
linv^oo ^i.2{t) — 0, but are otherwise arbitrary. One example is given by the stretched exponential law, or WeibuU 
distribution: 

$(t) = exp [ - [r(l + l/a)i^T]"l , < a < 1. (7) 

In lO, v = 1/(t) is a rate parameter having the meaning of inverse mean residence time (MRT) and r(a;) denotes the 
gamma-function. Moreover, the power law dependence, 

corresponds to the Pareto distribution. Both WeibuU and Pareto distributions typify the so-called fractal dependencies. 
In particular, such distributions have been detected for several different types of ion channels [Ullll. An interesting 
feature of the Pareto distribution is that for < 7 < 1 it displays a diverging variance, war(r) = 00, whereas the 
MRT (r) is finite. The closed time-intervals of a large conductance potassium ion channel studied in Ref. seems 
to obey ^ with 7 « 0.24. Other fractal-like distributions can be constructed from the expansion over exponentials 

oc 

*M ='^c^exp(-iy,T), = 1, (9) 

i—l i 

assuming some recurrence scaling relations among the rate constants {m;}, e.g., z^^+i = auj, ax\d among the expansion 
coefficients {ci}, e.g., q+i = ha, with some structural constants 0<a<l,0<5<l [l3,|23,|33- If the hierarchy 
of rate constants is obtained from a fundamental rate constant vq applying a recurrence scaling relation similar to 
one given above, the corresponding distribution can be characterized as a fractal in time. If the whole hierarchy is 
produced by a more complicated scaling law involving two, or more independent fundamental rate constants, the 
distribution is multi-fractal. The corresponding stochastic processes can be referred to as fractal renewal processes 
^2J. Such random processes presently attract renewed attention in physics and in mathematical biology |l2l |. 
The negative time-derivative 

^i.(r) = (10) 

yields the corresponding residence time distributions Next, let us assume that a number of alternations occurred 
before the starting time point and the considered process became homogeneous in time before the observation 
started at tg- Then, for such persistent, time-homogeneous process the RTDs of the first time interval tq = ti — to, 
must differ from V'i,2(r) [H [H US IH HI ■ Namely jll]. 



(n,2) 

where (T1.2) is given by Eq. JQ). The corresponding survival probability of the first residence time-interval reads 

im r'^i2{t)dt 

(n,2) 

Moreover, if to choose ^1.2(^0) ^ the stationary values, ^1,2(^0) — Pi*2; the considered persistent process is stationary. 
From Eq. it follows that the two-state renewal process (TSRP) can be stationary only if the two mean residence 
times, (ti) and (T2), are finite. A diverging mean residence time leads to anomalously slow diffusion (subdiffusion) in 
the multi-state case 113,133,131! such situations will not be addressed with this work. 

When a time-dependent input signal is switched on, the driven TSRP becomes a non-stationary process and the 
corresponding survival probabilities depend not only on the length of time intervals, but also on the initial time instant 
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t' of any considered residence time interval, i.e. ^i^2{t — t') ^i^2it,t'). The residence time distributions are then 
accordingly given by 



V'i,2(i,t') 



dt 



(13) 



The corresponding conditional survival probabilities can be defined as ^i^2{T\t') <&i,2(^' +T,t') (here the condition 
is different from that used in footnote [33 - in the absence of signal - notwithstanding the use of identical notations). 

The particular choice, $i.2(i,t') = exp ^ — J^* vi^2{'T)dT^ , leads to Eq. Q - the only choice which is consistent with 
the Markovian assumption 17]. In the non-stationary driven case, the distinction between ^'f\(t,t') and <i>i^2(i, i')j 
^'"^\{t,t') and 'ipi^2{t,t') is not necessary. Nevertheless, we keep formally this distinction in the following, because 
when the driving is being switched off, the process x{t) relaxes to its stationary state. This distinction becomes very 
important in order to construct the evolution operator for time-homogeneous initial preparations. 

B. Integral equations of non-stationary renewal theory 

Our immediate goal is to obtain the evolution equations for the considered stochastic process: we are looking for 
the forward evolution operator Tl{t\tQ) (or the matrix of conditional probabilities) connecting the probability vector 
p{t) = [pi{t),p2{t)]'^ at two different instants of time t and to; i.e.. 



Pit) = U{t\to)p{to) . 



(14) 



This evolution operator can be explicitly constructed by considering the contributions of all possible stochastic paths 
leading from p{to) to p{t). To start, let us separate these contributions as follows 



n(iito)-^n("H^|io), 



(15) 



n=0 



where the index n denotes the number of alternations that occurred during the stochastic evolution. The contribution 
with no alternations obviously reads. 



<^^°\t,to) 



<ff\t,to) 
Stochastic paths with a single alternation contribute as 



n(i)(t|to) = / dtiP(t,ti)F(")(ti,io), 

Jto 



where 



and 



F{t,to) 



$i(i,io) 
$2(i,to) 



4°^(t,to) 
^f\t,to) 



(16) 



(17) 



(18) 



(19) 



Next, the paths with two alternations contribute to Eq. (|15|1 as 



n(2)(t|to)= / dt2 / dtiP{t,t2)F{t2,ti)F^°^{ti,to), 
Jto Jto 



where 



F{t,to) 



Mt,ta) 
Mt,to) 



(20) 



(21) 



Contributions with higher n are constructed along the same line of reasoning. 
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This representation of the evolution operator Il(t\t') in terms of an infinite sum over the stochastic paths is exact, 
although not very useful in practice. The structure of the infinite scries in Eqs. p5|l - (|21|l implies, however, the 
following representation 

U{t\to)=n^°Ht\to)+ f dti-p{t,ti)G{ti,to), (22) 

J to 

where the unknown auxiliary matrix function G(t,to) satisfies the matrix integral equation 

G{t,to)=F^°\t,to)+ [ dhF{t,ti)G{h,to). (23) 

Jto 

The equivalence of Eqs. lfT5|) ~ (|^ and Eqs. can be readily checked by solving Ea. (|^ with the method of 

successive iterations. 

In components, Eq. H22|l reads 

nii(i|to) = <^^i \t,to) + f $i(t,ii)Gii(ti,io)dii , (24a) 



n22(t|io) = $r(i,io) + / ^2it,ti)G22itl,h)dti , 



(24b) 



ni2(i|to) = / <^i{t,h)Gi2{ti,to)dh , 
Jt„ 



(24c) 



n2i(t|io)=/ <^2{t,ti)G2i{ti,to)dti . (24d) 

It is worth to notice that the set of Eqs. (|24a|) - (|24d|) is not independent. The conservation of probability implies that 

Uii{t\to)+n2i{t\ta) = 1 , 

Tl22{t\to)+Tli2(t\to) - 1. (25) 

The consistency of Eqs. (|24af) - (|24dp with the conservation law, Eq. (|25|l . can be checked readily. The matrix integral 
eauation H23|l reads in components 

Gii(i,io)= / Mt,ti)G2iiti,to)dti , (26a) 

Jto 



G22{t,to)^ / Mt,ti)Gi2{ti,to)dti 

Jto 



(26b) 



Gi2(t,to) =4"^(i,io)+ / Mt,h)G22{ti,ta)dti , 

Jto 



(26c) 



nt 

G2i{tM ^ + MtM)Gii{tiM)dti . (26d) 

Jto 

From Eqs. H26a() - (|26d|l one can deduce independent scalar integral equations for each component of matrix function 
G(t,io)- Indeed, after substituting G2i(i,to) from Eq. Ij26d|) into Eq. (|26a|l the closed equation for Gii(i,to) follows 
as 

Gil (t, to) = ^f\tM + / a(i,ii)Gii(ti,to)dti . (27) 

Jto 
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In (EH), 

^[°\t,to)= I Mt,ti)i;^i\tuto)dh (28) 



and 



is a renewal density. Analogously, 



to 



a(t,to)=/ Mt,ti)Mh,to)dti (29) 



G22{t,to)^£.i°\t,to)+ f a{t,ti)G22{ti,to)dh , (30) 



to 



where 

(t,io) = / Mt,ti)lb^2\tl,to)dh , 
"'to 

6(^,^0)= / Mt,ti)Mti,to)dti . (31) 

Jto 

Moreover, for the ofF-diagonal elements of G{t,to) we find 

pt 

Gi2(t,io) = 4°^(t,to)+ / ^iit,h)Gi2{ti,to)dh , 



nt 

G2i{tM) - i'f\tM) + / 6(i,ti)G2i(ti,to)dii- (32) 

"'to 

Eqs. H27|l - (|31|l together with Eqs. H24a|) - (|24dp present a first main result of this work. This set of equations 
generalizes the integral equations of renewal theory obtained by Cox and others to the case of non-stationary 
renewal processes modulated by external signals. The solution of the evolution operator n(i|to) is thereby reduced 
to solve the set of independent scalar integral equations for Gij{t,to). This presents an essential simplification as 
compare to the case of an evaluation of infinite matrix integral series in Eqs. (|15|I - H21|I . 

C. Time-homogeneous case 

In the absence of a signal, all two-time quantities depend only on the time-difference t = t — to- In this case, the 
integral equations of renewal theory can be solved formally by use of the Laplace transform method and the evolution 
operator (i. e., its Laplace-transform) can be found explicitly. Let us denote the Laplace transform of any function 
F{t) below as F{s) :— exp{—sT)F{T)dT. Then, upon Laplace transforming Eqs. H24al) - (|32|l . using Eqs. 
(|12|1 and some well-known theorems of Laplace transform, one finds the explicit expression for the evolution operator 
n(s). It coincides with the known result in the literature [l3llT8ll2^ . reading 



fl{s) = - 

s 



^ _ G{sl Gjs) 

^ s{ti) s{t2) 

G(s) ^ _ G(s) 
s{ri) s{t2) 



(33) 



where 



(l-V'l(s)) (l-^2(s)) 

G{s) = ^ ^^^^ ^ (34) 

is an auxiliary function. 

The existence of finite mean residence times (ti 2) implies the following useful representation for the Laplace- 
transformed RTDs: 

V^i,2(s) := 1- (ti,2).s[1 + <?i,2(s)]. (35) 

In H35|l . (?i,2(s) are corresponding functions vanishing at s — > 0, i.e., 171,2(3) 0. Note that the functions 31,2(5) are 
not necessarily analytical. For example, g(s) ~ s'^ with some real- valued exponent, < 7 < 1, is allowed, for an 
example see below in Eq. (|94|l . Such non- analytical feature leads to diverging variance of RTDs. From the formal 
expression 1)33(1 a number of important results follows: 
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1. Stationary probabilities 



The vector of stationary probabilities p®* = [pf*,P2*]"^ '^^''^ be evaluated as p®* = linis^o (^511(5)^(0)^ With Eqs 
one readily obtains the result 



(ri) + (T2)' ^2 (ri) + (r2) 



i!. Relaxation function 

The generally non-exponential relaxation of {x{t)) — xipi{t) + X2P2{'t) to the stationary mean value Xst = a;ipf* + 
X2P2 is described by the relaxation function R{t), i.e. 

Pi,2(to + r) - p5 + [pi,2(<o) - Rir), (37) 
where i?(T) obeys the Laplace-transform 

m = - - f A + A) ^Gis) , (38) 

and G(s) is given by Eq. 1)34(1 . The validity of Eqs. 1)37(1 . 1(38(1 can be easily checked upon the use of Laplace- 
transformed Eq. ((14(1 and the result in Eqs. ((33(1 . ((34(1 along with the normalization condition pi(to) +^2(^0) = 1 
and Eq. ((36(1 . It should be emphasized here that that the relaxation function R(t) for the the considered persistent 
renewal process is unique, i.e. it does not depend on pi,2(io)- This corresponds to the situation where the random 
process x{t) has not been prepared at i = in a particular state Xi, or X2, but rather has almost relaxed to its 
stationary state. In other words, a number of alternations occurred before t — to and the probability pi,2(io) to 
measure the particular value xi^2 of x{t) at the instant of time to is close to its stationary value pf*2- This class 
of initial preparations, where the relaxation function does not depend on the actual initial probabilities, is termed 
the time-homogeneous preparation class. This preparation class |3,'il l3^ must be distinguished from strongly non- 
equilibrium initial preparations, where the system is prepared, for example, in a particular definite state, say in the 
state xi, with the probability one, pi(to) = 1- 



3. Stationary autocorrelation function and regression theorem 



Let us consider next the normalized autocorrelation function, i. e.. 



Hr) = hm (39) 



of the stationary fluctuations, 5x{t) — x{t) — Xst- In (|^ - 



(fe^)., = (Ax)^ ArrrA' (^o) 

((n) + {t2)Y 

is the mean squared amplitude of the stationary fluctuations and Aa: — X2 — xi is the fluctuation amplitude. With 

{5x{t + T)5x{t)) = {x{t + T)x{t)) — {x)lf., as t ^ 00, and 

lim (x{t + T)xit)) ^ y y X^XJU,J{T)pf (41) 
i=l,2j=l,2 



we obtain the same result as in Ref. j38| : i.e.. 



Hs)^--f^ + ^)\G{s). (42) 

S \{ti) {t2)J s2 



Upon comparison of 1(38(1 with ((42() we find the following regression theorem for these non-Markovian two-state 
processes; namely 

i?(r) = k{r) . (43) 
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The regression theorem (|43|) . which relates the decay of the relaxation function R{t) to the decay of stationary 
fluctuations klr), presents a cornerstone result for the derivation of phenomenological linear response theory for 
non-Markovian SR . 

Usually, the Laplace-transform (|42|l cannot be inverted analytically. If k{t) > for all times t, one can define the 
mean correlation time: 



k{t)dt = lim k{s 



(44) 



Assuming finite second moments of RDTs, (t^2) — Jo^ ''''^ "4^1,2 {T)dT we obtain from Eqs. (|42|l and 144|) the simple 
result 



where 



Rnmtm, 



(n)(r2) 
(n) + (r2) 



(45) 



(46) 



is the correlation time of the Markovian process possessing the same MRTs (ti_2} as the considered non-Markovian 
process. The coefficient 



Rnm — 2 (Ci + Q 



(47) 



presents a numerical quantifier of non-Markovian effects in terms of the coefficients of variation of the corresponding 
residence time distributions; i. e., 



Ci,: 



For example, for the stretched exponential Q the coefficient of variation emerges as 



C 



For the Pareto law distribution in © it reads 



C = 



I T{l + 2/a) 
r2(l + l/a) 



oo,7 < 1; 
^, 7>1 



(48) 



(49) 



(50) 



As a criterion for Markovian vs. non-Markovian behavior one can propose to test the coefficients of variation Ci^2 
of the experimentally determined RTDs V'i.2(i)- the strict Markovian case we have Ci = C2 = 1. Large deviations 
of any of the two coefficients of variation, Ci^2, from unity indicate the presence of strong non-Markovian memory 
effects. The proposed test-criterion appears experimentally to be more conveniently applied than the direct test of the 
Chapman-Kolmogorov-Smoluchowski equation |37| |. For example, in the fractal model of the ion channel gating by 
Liebovitch et al. the closed residence time distribution is fitted by Q with a w 0.2 [l^. This yields Cdosed ~ 15.84. 
Thus, assuming that the open residence times are exponentially distributed, i.e. Copen = 1, one obtains Rnm ~ 126. 
Furthermore, according to Ref. 21] BK ion channels display a closed residence time distribution following a Pareto 
law with 7 K, 0.24. In such a case, the memory effects should depict an infinite range since Tcorr = oo- In both cases, 
the observed two-state fluctuations do exhibit long-range temporal correlations. The gating dynamics is thus clearly 
non-Markovian within such a two-state description. 



4.. Power spectrum of fluctuations 
For the power spectrum of fluctuations; i.e., 

POO 

Sn{(^) = '2{Sx^)st k{t) cos{ujt)dt ^ 2{5x'^) stRe k{iuj) 
Jo 



(51) 
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xst(fo) 







FIG. 1: Relaxation of a perturbed persistent renewal process x{t). A constant force /o is applied long before and is released at 
t — to- The mean value {x{t)) relaxes from the constrained stationary value Xst{fo) to its true stationary value Xst- 



the use of Eqs. gDJ, (021) in (ED yields [H IM IM IM 111 



(n) + (t2) lu' 



:Re 



(52) 



It is evident that asymptotically, in the limit w — > oo, the power spectrum H52|) is Lorentzian in the case of time- 
continuous RTDs mila, 



Sn{uj) 



(Ti) + (r2)c.2' 



as — !■ oo . 



(53) 



This follows from the fact that limi^^oo '0i.2(*w) — and thus lim^^^oo G(ia;) = 1 39]. Practically this situation 
occurs for uj ^ ■ On the other hand, one can deduce from Eq. (|51|) that in the opposite limit for a; — > 0, 



(54) 



where {5x'^)st is the mean-squared amplitude of stationary fluctuations given by (|40|l and Tcorr is given in Eq. H45|l. 
A very interesting situation emerges for Tcorr — > 00, implying Sn{Q) — > 00. This occurs when at least one of the 
residence time distributions possesses a diverging variance, cf. Eqs. H45|l - H47(l . In such a case, for the low-frequency 
region uj < {ti^2)~^ the power spectrum drastically differs from the Lorentzian form. For example, for a symmetric 
TSRP with the survival probabilities given by the Pareto distributions ^ one can show (see also below) that for 
< 7 < 1, <S'Ar(tt') ~ l/uj^~''. For 7 — > this corresponds to celebrated 1// noise p3.l40|. 



III. PHENOMENOLOGICAL THEORY OF LINEAR RESPONSE 



It is possible to predict the linear response of the underlying stochastic process x{t) to the external driving f{t) by 
referring only to information on its stationary properties, i.e., without explicit knowledge of the concrete mechanism 
at work by which the process x( t) is pe rturbed by the external signal. The jjhenomenological theory of linear response 
for general stochastic processes and for thermal physical systems 42j provides a very useful and widely applied 

tool to answer this question. It is also the only method available if no further detailed knowledge of the microscopic 
dynamics is at hand for the observed two-state dynamics. This is the common experimental situation. The common 
linear response approximation 



(Sxit)) {x(t)) - Xst 



x{t - t')f{t')dt', 



(55) 



holds independently of the underlying stochastic dynamics jSg. In (|55|l . x(^) denotes the linear response function in 
the time domain. The universality of the relation (|55|l allows one to find the linear response function xit) using a 
properly designed form of the perturbation f{t). Within the phenomenological approach it can be obtained following 
an established procedure (i) First, apply a small static "force" /o, (ii) then, let the process x(t) relax to the 

constrained stationary state with mean value Xstifo), and finally (iii) suddenly remove the "force" at t = to, see Fig. 
□ 



Then, in accord with (|55|l the response function reads 



X{r) 



1 d 
fo dr 



{5x{to+T)), T>0, 



(56) 
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where {5x{to + t)) = xipi(to + r) + X2P2{to + ''') is determined by (|37() with the initial ^1.2(^0) taken as ^1,2(^0) = 
('''i,2(/o))/[('''i(/o)) + (■^2(/o))]- The hmit /o ^ is impUcitly assumed in Eq. (|^ . Expanding ^1^2(^0) to first order 
in /o we find with Ax = 2:2 — a;i 

(fe(to + r)) = ^^^[/32 - /3i]i?(r)/o + o(/o), (57) 



where 

. rfln(ri.2(/o)) 

Pl,2 •— -Ti 



/o=o . (58) 



Note that in the derivation of this result it is tacitly assumed that the initial constrained stationary populations 
Pi, 2(^0) at t = to belongs to the class of time-homogeneous initial preparations ^6J for the process x{t) in the absence 
of applied force. This seems a natural and intuitively clear assumption in view of the facts that the limit /o ^ has to 
be taken in H56|l at the very end of calculation, and the considered process is persistent. Nevertheless, this commonly 
accepted assumption is a hidden hypothesis which, strictly speaking, cannot be proven within the phenomenological 
approach. 

Upon combining 157|) with the regression theorem (|43ll we obtain from H56|l . after taking the limit /o 0, the 
fluctuation theorem ,22j 

X{r) = -[/32 - + r)5x{t)),, , (59) 

wherein 0{t) denotes here the unit step function. The non-Markovian fluctuation theorem (|59|l presents a prominent 
result 28]; in particular, it does not assume thermal equilibrium |3^ . In the frequency domain it reads 

X(w) = ^1 + «wfc(-iw)J , (60) 

where x{'^) = x{^)^^^^dt denotes the linear response function in the frequency domain, and k{s) is given by Eq. 
Substitution of Eqs. lO and ^ in (jHOJ) yields 

{(32 ~ /3i)Aa: i . 
Xi'^) = / \ , / N -G{-iuj), (61) 

where G(s) is given in 134(1 . The expression H61|l together with Eq. ((34|l connects the linear response function x{Ld) 
with the Laplace-transformed residence time distributions ?/'i,2(^i^), i-e., with the characteristic functions of the RTDs. 
If, in addition, the mean residence times obey the thermal detailed balance relation 

{r,[h)) _^^f^^n^Ml)^ (62) 



(r2(/o)) V keT 

where e(T) is the free-energy difference between two metastable states, we recover for the fluctuation theorem in l|59|) 
the form that characterizes classical equilibrium dynamics [3^E3.E3; i-^-, 



Eq. (|nU then yields 



x{r)^-^-^^{5x{T)5xm,,. (63) 



= ^^""l \ -G{~^u;) . (64) 
ksT (ti) + (T2) uj 

For example, this result is valid for an Arrhenius-like dependence of (ti_2) on temperature T and force /o; i.e. 

/ ffw A /AC/1,2 TAa;i,2/o\ ^^.x 
(n,2(/o)> = ^i,2exp(^ — j, (65) 

where AC/1^2 are the heights of activation barriers, Axi = zAx, Ax2 = (1 — z)Ax with Aa; = X2 — a^i, < z < 1. 
Eq. presents a key result because it provides a link between the phenomenological theory of linear response 

theory and the actual physical processes which are in thermal equilibrium and do exhibit long-range time correlations. 
Let us assume, for example, the following situation: The observed two-state process results from thermally activated 
transitions in a complex potential energy landscape U{x) possessing two domains of attraction (i.e., two metastable 
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states) separated by distance Aa; along the direction of the reaction coordinate x which describes transitions between 
the metastable states. Next, let us assume that the coupling of the external force f{t) to the dynamics has the potential 
energy form Uint = —xf{t). Then, the classical equilibrium fluctuation theorem (jfiHfl follows from first principles [i^ . 
or, likewise, from a mesoscopic starting point in terms of the generalized master equation for the thermal equilibrium 
dynamics |43| : in other words, it is exact. The non-exponential features of the RTDs in the described situation 
stems from the motions "perpendicular" to the above reaction coordinate x. In such a case, the thermodynamic 
relations like Eq. Ht)2|) are compatible with non-Markovian kinetics. This is the case where the phenomenological 
theory of linear response in non-Markovian systems has a firm foundation. The readers should be warned, however, 
that the phenomenological theory is not universally valid for nonequilibrium physical systems. Nevertheless, below 
we explicitly define an universality class of such systems (which are beyond the thermal equilibrium class) where its 
validity can be proven on a more general basis. 



IV. ASYMPTOTIC RESPONSE THEORY BASED ON DRIVEN RENEWAL EQUATIONS 

Starting from the driven renewal equations H24a|l - (|32|l one can develop the theory of the linear and the nonlinear 
response which possesses a broader range of validity as compared to the above phenomenological theory. For a periodic 
signal (switched on in the infinite past) like in jSJ, the conditional survival probabilities $i,2('''|^) ■— ^i.2{t + T,t) 
acquire (at asymptotic times t 3> to) the time periodicity in t of the driving signal and therefore can be expanded 
into the Fourier series, i.e., 

^iAr\t)= $l:2^(T)expHnm], 'pijHr) ^ [^^^^(r)] . (66) 

n— — OG 

Similar expansions hold also for the conditional residence time distributions V'i,2(''"|t) with the corresponding expansion 
coefficients ''Pi^2(j) ~ ^"^^ 1^2 (''')■ Note that ^flir) and ip^lir) in this section denote the Fourier expansion 
coefficients with n = 0. These quantities are clearly not related to the survival functions (|12|l and RTDs (|ll|l of the 
first time interval. We hope that such use of notations will not confuse the readers. The corresponding Laplace- 
transformed quantities of the r-dependent Fourier coefficients '0j"2(*) ^'^'^ ^i^2i^) related by 

^S(s) = '5„,o"sl>l"2^(s). (67) 

Our goal is to evaluate the asymptotic behavior of the populations (t) and of the mean value (x'-^'^t)). To do 
so, one needs to determine the asymptotic evolution operator n("*^(i) := limtQ_,_oo n(t|to)- Obviously, Il^^^'^^t) — 
nil'^t) and Il''22'\t) = ni^T\t)- Moreover, p['''\t) = n^i°^^(t), pi'"'\t) = n^^'Ht). Next, let us define the auxihary 
quantity G^'''\t) as G^'^'^t) := limto^_oo G{t,to). Then, Eqs. ifj^ . lf?7j) in the limit to -00 yield 



J — 00 



(68) 



where G^^^'' (t) is solution of the integral equation: 



G[r\t) - / ut,h)G[r\h)dh, (69) 

J — 00 

with the renewal density given in Eq. (|29|) . The equation determining p^2^^\t) hkewise reads 

pi'^'Ht) = f ^2{t,h)Gif{h)dh, (70) 
where 0^22'' (^) ^h*^ solution of integral equation 

Gi7\t) = f Ut,h)Gir\h)dh, (71) 



with ^2(^,^1) given in Eq. H31|) . Note that the conditional renewal densities ^1,2(^1^) ■— £,i,2{t + T,t) also acquire 
a time-periodicity in t and can be represented in the form like H66|l with the corresponding expansion coefficients 
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Ci"2 (''')■ One can show that the corresponding Laplace-transformed quantities ^^"2 (*) related with the quantities 
''Pi^'ii^) follows: 

00 

m— — 00 
00 

ii"\s) ^ i>['-\s)4,i"-'-\s + ^mn). (72) 

m — — 00 

For periodic driving /(i), both i^) ^11^22(0 "lust be periodic functions of time 01 and can be expanded into 
Fourier series: 

p^r\t) = J2 Puh-'"'''^ pt2^ = iPi}r (73) 

A;— — 00 

and 

00 

G'r^Ut) = E ai'le-^"'''^ ait' = [9^ilr^ (74) 

fc— — oo 

respectively. 

Using Eqs. and the expansion (|73|l one can show that the coefficient in Eq. H73|l determines the linear 

response Junction x{^) in the frequency domain as 

m) - -'-^P^' (75) 

Jo 

in the limit /o 0. Moreover, from the normalization condition p^'^*'' (t) + p^2^' {^) = 1 it follows that 

Pi' + P2' - 1, P^r' = -Pt' for n / 0. (76) 

Upon substituting Eqs. H73|l and 174() and the expansions like Eq. IjGGI) into Eqs. (|68|I - H71|I . performing the time inte- 
gration and comparing the cocfhcicnts of the Fourier expansions on the left and right hand sides of the corresponding 
equations we finally end up with: 

00 

Pf' = E ^i"'(-*fc^^)fff'"\ (77) 

n— — 00 

00 00 

- E E i^'^\-^km^r"'\^^[k-mm9f-''\ (78) 



n— — 00 m— — 00 



9i 
and 

00 

P^2' = E ^t\-^m9t"\ (79) 

71— — OG 

oc 00 
n— — cxD m— — 00 

The relations l|77|) - (|80|l also serve as the basis for a response theory without restriction on the linear response approx- 
imation. In order to apply these equations, one has to specify the expansion coefficients in Eq. ]<M\i . i.e. to concretize 
the way how the external signal f(t) enters the conditional residence time distributions 'ipi^2{T\t), or, equivalently, 
the conditional survival probabilities $i,2(t|0 to the required order in the signal amplitude /q. It is worth to notice 
that the solutions of Eqs. (|78|) . H8U|) are defined up to some arbitrary constants which can be fixed at the end of 
calculations by applying the normalization relations in (|76|l . 

In the linear response approximation, $x*'2('^) ~ '^'i.2('S), i-e., ^^fli^) coincide with the unperturbed survival proba- 
bilities $1,2(3). Moreover, ^'"ilis) oc /q. All the higher order terms ^i^^^^'is) can be neglected, being of higher order 
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proportional to Jq, n > 2. The same holds true for "01 2 (*)- After some cumbersome algebra, one finds from Eqs. 
(|77|l - (|80ll an expression for p\^\ which then by use of relation H75|l yields 

The result in Eq. (|81|) presents a second cornerstone result of this work. Note that this general result depends on the 
quantities ip[^l{s) oc /o which do not follow directly from the characteristic functions of stationary RTDs, i.e., V'i.2(s), 
but their knowledge requires one to specify a microscopic model. Generally, Eq. (|81|) is not mathematically reducible 
to the result H61(l of the phenomenological theory. A question arises whether such a reduction is possible in practice 
and the phenomenological theory of linear response can be put on a more firm ground beyond the time-homogeneous 
preparation class result in Eq. H(il|l of which the thermal equilibrium result in Eq. I|64|l is a special case. Below we 
describe a rather broad class of relevant systems. 



A. Models with form-invariant RTDs 



Let us assume that the survival probability and the corresponding RTD can be parameterized by a single frequency 
parameter v with has the meaning of an inverse mean residence time, i.e., v = (t)^^. Furthermore, we assume that 
a weak signal f{t) causes v to became time-dependent, i. e.. 



(82) 



with [3 ^ l//o (the subscripts 1,2 are suppressed). Moreover, the survival probabilities become modified applying 



the following rule: vt Jj*^^ v{t')dt' . More generally, let us consider arbitrary survival probabilities of the form 



generalized to the time-inhomogeneous case in the following way 

(•t + T 



$(r) ->$(r|t) =^c,exp(^- j v,{t')dt'y ^ 1. 



(83) 



In 183|l . we assume that (to leading order) neither the expansion coefficients Ci, nor the ratios between any of Vi{t) 
and Vj{t) are modified by the applied signal f{t), i.e.. 



(84) 



with aij being some structural constants. This covers fractal (although not multi-fractal) time distributions. Put 
differently, the scaling law by which the whole hierarchy of rate constants is produced from a single rate constant is 
invariant of the applied signal. If the mean residence time (r) — Ci/vi exists, one can always set v — (t)^^ as the 
relevant rate constant in the absence of driving. This rate will acquire an explicit time-dependence like in (|82|l when 
the signal is switched on. Given our assumptions, all the time-dependent rates Vi{t) in H83(l will be proportional to 
the rate v{t) in Eq. (|82|l . Then, in the lowest first order in /3/o, we find 



$(T|t) = $(t) -f /3V(t) 
From ()85|l we obtain upon observing Eq. (O 



t + T 



f{t')dt' 



1 ./9l,2/o , / X r I -n, \ 1 

2«— ^V'i,2(T)[exp(-iS2T) - 1 



and 



l, /?l,2/0 

2* n 



-01,2(5) - 01, 2(5 + iO) 



(85) 



(86) 



(87) 



Observing Eq. Ht)7|l by taking into account 0i.2(O) = 1 in Eq. (|87|) thus yields 



:/3i,2/o[l-0i,2(-«f^)] 



Substituting H88(l into H81|l we recover the result of the phenomenological theory in Eq. H61|l . In conclusion, for 
the considered class of models the non-equilibrium fluctuation theorem Ijtilf) is well justifled. This model class can 
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therefore be reconciled with the assumption of time-homogeneous initial preparations used in the phenomenological 
theory of linear response (see Sec. III). This assumption is naturally not always justified a priori. It rather delimits an 
important and rather broad class of corresponding physical systems. Nevertheless, the equilibrium fluctuation theorem 
(|63|1 presents a fundamental relation which must be obeyed for all thermal equilibrium systems. This imposes a salient 
restriction on mesoscopic models leading to the observed equilibrium non-Markovian dynamics. In particular, if one 
knows that the considered system is in the thermal equilibrium, one must use the rigorous relation (|64|l . rather than 
(|81|l for the calculation of the linear response. This constitutes the essence of the phenomenological theory of non- 
Markovian stochastic resonance developed in Ref. 28]. For other systems, e.g., for those modelling neuronal dynamics 
(which are far away from thermal equilibrium) the use of Eq. H81|l is preferred. In order to apply Eq. (|81|l . however, 
one must also specify the underlying non-equilibrium microscopic dynamics in the presence of a time-periodic stimulus. 
This means that the time-inhomogeneous conditional RTDs i/ji 2(T|t) must be measured, or modelled (to the linear 
order) in the driving signal strength. We next present a detailed study of non-Markovian Stochastic Resonance in 
thermal equilibrium systems that do exhibit prominent temporal long-range time correlations 28] . 



STOCHASTIC RESONANCE 



In the presence of applied periodic signal (jSJ, the spectral power amplification (SPA) 0,0], r]{Q) 
by use of the fluctuation theorem in upon combining , , , as follows 



{kBT? cos\v^[e{T)l{2kBT)] 



\x{n)\^ reads 



(89) 



In (jHHl), i^iT) = (ti}"1 -I- (ra)"^ denotes the sum of effective rates. The quantity e{T) = AU2 - MJi + TAS denotes 
the free-energy difference between the metastable states which includes the entropy difference AS :— S2 — Si = 
kBin{A2/Ai). In the Markovian case we obtain G(s) — s/{s + v) and (|89|l equals the known result, see in [^.lil]. 
The signal-to-noise ratio (SNR) is given within linear response approximation by 



SNR(r2,r) := 



SNm ' 

where Sn{uj), Eq. (|51() . is the spectral power of stationary fluctuations 2]. By use of H89|l . we obtain 

7Tf^{Ax/2)^ v{T) 



SNR(r2, T) ^ 



2kBT} 



where the term 



N{n) = 



Re[G{in)] 



(90) 



(91) 



(92) 



denotes a frequency- and temperature-dependent non-Markovian correction. For arbitrary continuous "01.2 (t) the 
function A^(Sl) approaches unity for high-frequency signals, SI ^ (■'■1,2) Thus, Eq. H91|) reduces in this limit to 
the known Markovian result 0, i.e. the Markovian limit of SNR is assumed asymptotically in the high frequency 
regime. More interesting, however, is the result for small frequency driving. In the zero-frequency limit we find 
A^(0) = 1/Rmm with Rmm given in Eq. (|47|l . With Rmm = cxo as it is the case for the Pareto distribution (jSJ with 
< 7 < 1, N{{)) = 0; i.e. SNR(17 = 0,T) = as well. Consequently, ultra-slow signals are difficult to detect within 
the SNR-measure in a strongly non-Markovian situation. 



A. Symmetric SR 

As a first example, we address non-Markovian SR in a symmetric system with the survival probabilities <I'i,2(''') 
described by the identical power laws lO with v = (t)^^ determined from Eq. Q with f{t) = 0, = and 
Ai7i 2 = AU. In this case, the Laplace-transformed RTDs read 

i^{s) = 1- (7(r)s)''^'exp(7(r>5)r(-7,7(T)s), (93) 

where T{x,y) is the incomplete gamma-function |46j. For < 7 < 1, the distribution (|93() has a diverging variance; 
its small-s expansion reads 

i^{s) « 1 - {t)s + 7^r(l - -f)[{T)s]^+'>. (94) 
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FIG. 2: The spectral amplification of the signal (in arbitrary units) is depicted versus the thermal noise intensity D — ksT at 
different driving frequencies Q.: (a) non-Markovian symmetric system and (b) its Markovian counterpart. In the non-Markovian 
case, both RTDs follow a Pareto law with 7 = 0.2. D is scaled in units of the barrier height A[/; n is scaled in units of vq. 



Using (|94() and H34() in H52|l we obtain for the low-frequency part of the power spectrum 

Sm{oj) « i(Aa;)^r(l -7)sin(7r7/2)^;^^2(^. (95) 

To obtain the spectral amplification H89|l and the SNR H91I) numerically one has to use Eq. (|93|1 in Eq. (|34|l . For 
7 > 1, the power spectrum of this process mimics a conventional Lorentzian. Moreover, for 7 ^ 1, C « 1, cf. H5U|I . 
Thus, one can expect that for large 7 the considered situation does not differ much from the Markovian case, at least 
qualitatively. Indeed, for very large 7 ^ 100 the discrepancy with the Markovian case in the SNR behavior versus 
noise intensity D = ksT is not detectable. The well known bell-shaped Stochastic Resonance behavior is reproduced 
with the maximum at D = AU/2. Nevertheless, in the behavior of r/(il,T) some discrepancy still remains detectable 
even for such large 7 (not shown). 

Next, the case with < 7 < 1 is of major interest as it is qualitatively very distinct from the Markovian Stochastic 
Resonance, see Fig. 2. The reason is that the mean correlation time Tcorr in Eq. (|44(l becomes formally infinite and 
the power spectrum exhibits a typical l//"-characteristics, with a = 1 — 7, cf. Eq. (|95|l . Nevertheless, an important 
time scale of the stochastic dynamics does still exist: It is defined by the mean time of stochastic turnovers between 
the metastable states, to{D) = 2(t). Invoking the reasoning of a stochastic synchronization of stochastic resonance 
[4^ one can expect Stochastic Resonance to occur when the time scale of stochastic turnovers tq{D) matches the 
period of external driving T = 27r/r2, i.e., to{D) ~ T. Indeed, Fig. 2 (a) unambiguously demonstrates the Stochastic 
Resonance phenomenon for a non-Markovian system with 7 = 0.2. Thus, the interpretation of SR as the phenomenon 
caused by stochastic synchronization between the time-scales of the random, temperature driven transitions and the 
external periodic modulations 2, 44] can be extended even onto this extreme non-Markovian case (with diverging 
mean correlation time, Tcorr — 00). Note, however, that the maximal value of the spectral amplification of signal 
is strongly suppressed in the present case by the factor of about 20 as compare with the corresponding Markovian 
counterpart possessing the same (t), see Fig. 2(b). 

In contrast to the overall simpler behavior of the spectral amplification measure the SNR displays prime new 
features, cf. Fig. 3 (a) and Fig. 3 (b). First, the SNR becomes frequency-dependent. In the limit ^ 0, we obtain 
for the form-factor N{n) in Eq. 

2.in(4;2h'r(l-7) - 

In this limit, the signal-to- noise ratio can be approximated as 



(97) 
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FIG. 3: Signal-to-noise ratio (in arbitrary units) versus thermal noise intensity D — ksT at different driving frequencies SI; (a) 
non-Markovian symmetric system and (b) its Markovian counterpart. In the non-Markovian case, both RTDs follow a Pareto 
law with 7 = 0.2. D is scaled in the units of barrier height AU, fl is measured in units of vq- 



This SNR expression (|97() displays several nontrivial features: (i) the Stochastic Resonance peak occurs at the noise 
strength Dnm{^ ^ 0) = 7AC//2 as compare to the Markovian case, where Dm = AC//2. (ii) The SNR displays a 
nontrivial, power law dependence on the angular driving frequency SNR(f^) - n'^-^ . Moreover, with the increasing 
angular frequency of signal the signal-to- noise ratio, SNR(r2), should approach its frequency independent Markovian 
limit. Thus, the resonance value Dnm{^) becomes frequency-dependent for an intermediate range of frequencies and 
approaches the Markovian value Dm in the limit of high frequencies. This profound frequency-dependence of non- 
Markovian Stochastic Resonance is very distinct from its Markovian counterpart, compare Fig. 3(a) with Fig. 3(b). 



B. SR in ion channels with fractal kinetics 



Our second example pertains to the non-Markovian SR in an asymmetric system. An especially interesting case 
emerges when one of the RTDs is exponential, while the one presents a power law with a giant (divergent) dispersion. 
Interestingly enough, such a case apparently is realized in nature for the gating dynamics of the locust BK channel 
plf. Indeed, this and some other ion channels exhibit a fractal gating kinetics together with the 1//" noise power 
spectrum of fluctuations 21, 22i> 2di> In the context of gating dynamics, x{t) corresponds to the conductance 
fluctuations and the forcing f{t) is proportional to the time- varying transmembrane voltage. For a locust BK channel 
the measured unperturbed closed time statistics iPi{t) can be approximated by a Pareto law ||SJ) with 7 « 0.24 and 
(ti) = 0.84 ms |2lj. The open time RTD assumes an exponential form with (T2) = 0.79 ms 21]. 

Unfortunately, neither the voltage, nor the temperature dependence of the mean residence times are experimentally 
available. For this reason, we employ here the common Arrhenius dependence in (|65l) with some characteristic 
parameters. Namely, because the temperature dependence of open-to-closed transitions is typically strong 0, we 
assume a high activation barrier; i.e. AU2 = 100 kJ/mol (^ 40 kBTroom)- The closed-to-open transitions are assumed 
to be weakly temperature-dependent with AUi — 10 kJ/mol. Because (ri) ^ (T2) at room temperature Troom, the 
difference between At/i and /S.U2 is compensated by an entropy difference A5 ^ —36 fc^. The physical reasoning 
is that the closed time statistics exhibits a power law; i.e. the conformations in the closed state form a self-similar 
hierarchy and are largely degenerate j4^]. This in turn implies a larger entropy as compared to the open state. 

The normalized autocorrelation function k{i) and the power spectrum SNi^^) of the current fluctuations are of 
prime interest. In the considered case, the auxiliary function (|34|) simplifles to 

G(.) = ^^lMi^i^, (98) 

where "01(5) is given by Eq. (|93|l with (r) — The Laplace-transform of k{t) can in the limit s — > be approximated 
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r [msec] uj [sec ^] 



FIG. 4: (a) The normalized autocorrelation function of current fluctuations, see Eq. I39II . and (b) the corresponding power 
spectrum for the studied model of locust BK channel. The amplitude of current fluctuations is taken to be 10 pA. The broken 
line in (a) corresponds to the long-time asymptotic, Eq. IIUUII . being in agreement with the numerical result (full line) in 
long-time limit. 



as 

From (|99ll the long-time {t oo) behavior of the autocorrelation function follows immediately by virtue of a Tauberian 
theorem ,3Q|], namely 

fcW^P2'(:^) \ (100) 

where p|* is the channel's stationary opening probability. The result in H1UU|) describes a power law decay with an 
exponent 7 = 0.24. In Fig. 4(a), this analytical result is compared with the numerical inversion of k{s) with G{s) in 
(|98|l . obtained due to the Stehfest algorithm This figure shows that the long-time asymptotical behavior of A;(t) 
indeed obeys the power law in (|100() for r > 10 sec. However, for smaller r < 10 sec some kind of transient behavior 
occurs which cannot be characterized by a simple power law. Nevertheless, the slow decay of correlations is clearly 
non-exponential. 

For uj 3> (ti_2)~^ the power spectrum of fluctuations is expected to approach a Lorentzian tail, S{u!) ^ lo^"^ . Indeed, 
this behavior starts in Fig. 4(b) for oj > 500 sec~^. The nontrivial frequency dependence emerges for the sufficiently 
small cu ^ {ti.2)~^- In this case we obtain from l|99|l 

Sic ^ 0) « 2(Ax)2 Jlll!M!^r(l - 7)sin(7r7/2) ^ (101) 

Thus, for 7 — 0.24 we have S{uj ^ 0) cx l/w" with a — 1 — 7 = 0.76. This typical 1//" noise behavior is depicted in 
Fig. 4(b). We should remark, however, that the experiment |23| gives a slightly different value of a « 1. The reasons 
of this discrepancy are presently not clear. One possibility is that the durations of the subsequent open and closed 
time intervals are yet mutually correlated, contrary to the assumptions made in the present model. If this is the case 
indeed, the studied model should be generalized further to account for such correlations. 

The spectral power amplification versus the temperature is depicted for various angular driving frequencies in Fig. 
5(a). The panel in Fig. 5(b) corresponds to an overall Markovian modeling with an exponential iPi{t) possessing the 
same mean residence time (ti). We observe a series of striking non-Markovian features in Fig. 5: (i) A characteristic 
SR-maximum occurs in the physiological range of varying temperatures. This maximum is caused by entropic effects 
which have been not addressed before in the theory of stochastic resonance. Because of the fact that the free-energy 
bias e(T) is temperature-dependent, due to a large entropic asymmetry between states. Stochastic Resonance in the 
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FIG. 5: (a) The spectral power amplification ri{n), Eq. 1891 . (in arbitrary units) vs. temperature (in °C) for the BK ion 
channel gating scenario (see text) and (b) its comparison with a corresponding Markov modeling. 
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FIG. 6: The signal-to-noise ratio (SNR in arbitrary units) versus temperature (in °C) for the studied model of Stochastic 
Resonance in a locust BK channel. The upper curve depicts the Markovian limit attained for large angular driving frequencies 
of the signal. 



spectral amplification occurs in a temperature regime where the populations of both states become approximately 
equal, e{T) w 0. Note that this effect occurs also in the Markovian case, cf. Fig. 5(a) and Fig. 5(b). Therefore, 
it is not caused by non-Markovian effects, (ii) Due to a intrinsic asymmetry the (angular) frequency dependence of 
the spectral amplification 77(0, T) for the Markov modeling is rather weak for small frequencies V, <^ ('''1^2)""'^ Q- 
contrast, the non-Markovian SR exhibits a distinct low-frequency dependence, thereby frequency-resolving the three 
overlapping lines in Fig. 5b. This feature constitutes an authentic non-Markovian effect, (iii) The evaluation of 
the SNR yields - in clear contrast to the frequency-independent Markov modelling - a profound, very strong non- 
Markovian SR frequency suppression of SNR towards smaller frequencies: The SNR-maximum for the top line in 
Fig. 5(a) is suppressed by two orders of magnitude as compared to the Markov case, cf. Fig. 6. As a consequence, 
for a strong non-Markovian situation it is preferable to use low-to-moderate frequency inputs in order to monitor 
non-Markovian Stochastic Resonance with SNR. 
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VI. SUMMARY AND CONCLUSIONS 



In the present work we have put forward a general theory of Stochastic Resonance for two state non-Markovian 
systems. The theory is based on time-inhomogeneous integral renewal equations governing the evolution of conditional 
probabilities in the presence of driving signal. These equations for driven renewal processes generalize earlier result 
by Cox [3 and others for stationary renewal processes to include the signal influence on the residence time 
distributions. Based on these new equations we presented a general outline of the theory of the linear and the 
asymptotic nonlinear response to the sinusoidal signal. In particular, we obtained a general expression for the linear 
response function Eq. (|81|l . which can be used for a variety of applications. The expression in (|81|l presents 

major result of this paper. We note, however, that the explicit use of Eq. 181(1 requires one to specify explicitly 
the way in which the periodic signal modulates the asymptotic, non-equilibrium residence time distributions. For a 
class of non-equilibrium fractal distributions where the signal enters the RDTs through a single frequency parameter 
having the meaning of the inverse mean residence time, it has been shown that Eq. H81|l reduces to the result H61|l 
of the phenomenological theory of linear response developed previously in Ref. p^ . Moreover, if the mean residence 
times obey the thermal detailed balance condition (|62|l . the expression (jtill) reduces further to Eq. 1)64(1 which can 
be obtained independently from the classical fluctuation-dissipation theorem ((63(1 by use of the expression in 1(421) for 
the autocorrelation function of the considered non-Markovian stochastic process. Even though the microscopic (or 
mesoscopic) details of the thermal equilibrium dynamics leading to the observed two state non-Markovian fluctuations 
are generally not known, the linear response function is determined uniquely by the characteristic functions of the 
residence time-distributions V'i,2(s) via Eqs. and For such equilibrium non-Markovian fluctuations, the 

knowledge of the equilibrium RTDs allows one to determine the linear response of the considered physical system to 
weak signals. This is the essence of the phenomenological theory of non-Markovian stochastic resonance put forward 
in Ref. [S^. For such equilibrium systems, the general expressions for the spectral power amplification, Eq. 1(89(1 . 
and for the signal-to-noise ratio (SNR), ((91(1 . are available. We applied these general expressions to study the main 
features of Stochastic Resonance in several non-Markovian systems exhibiting long-range temporal correlations along 
with 1//" power spectra of fluctuations. 

In particular, for a symmetric non-Markovian system with a power law distributed residence time intervals the 
occurrence of Stochastic Resonance has been demonstrated to comply with a stochastic frequency synchronization 
similar to the Markovian case However, both the SPA- measure and the SNR- measure become strongly suppressed 
due to strong non-Markovian effects. The most striking feature of the non-Markovian SR is a distinct frequency 
dependence of the SNR measure. In particular, the SNR becomes immensely suppressed for low frequency signals. 
Thus, the use of signals with an intermediate frequency range matching the mean time of the stochastic escapes 
between states yields most distinct non-Markovian SR-feature. 

For asymmetric non-Markovian fluctuations pertinent to fractal gating dynamics of the locust BK ion channel several 
new features have been revealed, (i) The expected diminishment of the SPA-measure relative to the Markovian case 
does not occur. This can be attributed to the fact that one of the RDTs in the considered case is strictly exponential 
similar the Markovian case, (ii) For asymmetric Markovian systems the SPA-measure ceases to be frequency-dependent 
for small adiabatic frequencies. The non-Markovian effects, however, introduce at low driving frequencies a distinct 
dependence, both for the SPA and the SNR. This latter phenomenon can be used to detect and establish a strong 
non-Markovian behavior in practice. 

Our novel non-Markovian theory of Stochastic Resonance possesses a whole range of applications and we hope that 
it will be used by the practitioners in their further research work on Stochastic Resonance. Especially, we hope that 
our theory will guide experimentalists to find the proper and most interesting parameter regimes and to reveal the 
Stochastic Resonance effect on the level of single biomolecules. 
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